function [tau1, tau2, offset, a1, a2] = multitrace_fit()
% multitrace_fit.m - fit multiple traces with shared tau's, double exponential
% 10/25/2001 Paul B. Manis, Ph.D.
% pmanis@med.unc.edu

% access current data.
d = datac('getdfile');
y_data = datac('geti');
x_data = make_time(d);
v = datac('getv');
v = v(:,[1000:1500]);
vm = mean(v');

delay = floor(5 * d.rate(1)*d.channels(1));

% reduce number of points.
pts = [2250:4900];
% 2250+unique(floor(record_parse('1;2700/250l')));

[tau1, tau2, offset, a1, a2] = mtf(x_data(:,pts), y_data(:,pts), vm, 500);

% plot voltage-dependence of amplitudes.
h = findobj('tag', 'multiexp');
if(isempty(h))
   figure('tag', 'multiexp');
else
   figure(h)
end;

subplot('position', [0.1 0.1 0.8 0.35]);

plot(vm, a1, 'rx');
hold on;
plot(vm, a2, 'bs');

return;

function [tau1, tau2, offset, a1, a2] = mtf(x_data, y_data, vm, maxiter)
% routine thata actually does the work.
% fit multiple traces simultaneously
% with same taus (two exponential fit...).



u = size(y_data);

x_data = x_data - min(min(x_data)); % offset to 0 time
alpha = 0; beta = 0;

% first, fit some traces individually to get initial parameter estimates

w = find(vm <= -80);  % limit to traces hyp beyond -80 mV.

for I = w
SIG = [];
VP = [1 1 1 1 1]; % floating TAUS.
LB = [-100 0 1 0 1];
UB = [5000 15000 50 15000 500];
   c = max(y_data(I,10:end));
   [FP, CHISQ, NITER, fx] = MRQFIT('exponential', [0 c*0.7 10 c*0.3 300], x_data(I,:), y_data(I,:), ...
      SIG, VP, LB, UB, 100);
   %   s_model = 51; % model in fit_func.c : double exponential.
   %   s_initpar = [0 c*0.7 10 c*0.3 300];
   %   s_order = length(s_initpar);
   %   s_pmask = [1 1 1 1 1];
   %   s_lbound = [-1000 0 2 0 75];
   %   s_ubound = [1000 10000 75 10000 500];
   %   [c, s_lam] = curve_fitting(x_data(I(1),10:end)', y_data(I(1),10:end)','simplex','cubic', s_model, s_order, s_initpar,...
   %      s_pmask, s_lbound, s_ubound, alpha, beta, 100 );
   
   % compute the fit trace  
   %   fy(I,:) = -fit_func(s_lam,x_data(I(1),:),0*y_data(I(1),:),s_model,alpha,0);
   fy(I,:) = fx;
   
   % store the results
   tau1e(I) = FP(3); % s_lam(3);
   tau2e(I) = FP(5); % s_lam(5);
   amp1e(I) = FP(2); % s_lam(2);
   amp2e(I) = FP(4); % s_lam(4);
   os(I) = FP(1); % s_lam(1);
end;


% now fit all of the traces, using the mean taus from the first fits.
tau1x=find(tau1e ~= 0);
tau2x=find(tau2e ~= 0);
tau1e = tau1e(tau1x)
tau2e = tau2e(tau2x)

% taus
initpar(1) = mean(tau1e); % s_lam(3);
initpar(2) = mean(tau2e); % s_lam(5);

for i = 1:u(1)
   c = max(y_data(I,10:end));
SIG = [];
VP = [1 1 0 1 0]; % FIX TAUS.
LB = [-100 0 1 0 1];
UB = [5000 15000 50 15000 500];
   [FP, CHISQ, NITER, fx] = MRQFIT('exponential', [0 c*0.7 mean(tau1e) c*0.3 mean(tau2e)], x_data(i,:), y_data(i,:), ...
      SIG, VP, LB, UB, 100);
   z(i,:) = fx; % save data   
   tau1(i) = FP(3); % s_lam(3);
   tau2(i) = FP(5); % s_lam(5);
   a1(i) = FP(2); % s_lam(2);
   a2(i) = FP(4); % s_lam(4);
   offset(i) = FP(1); % s_lam(1);
end;

if(exist('oldcode'))
   %----------------------------
   % now fit all of the traces. Tau's are fixed from the determination above and are
   % shared among traces. We allow the amplitudes and the offset values to vary for each trace.
   %
   % seed initpar with lam's from above. Note order is different for 51 vs 74 (fit_func.m)
   tau1x=find(tau1e ~= 0);
   tau2x=find(tau2e ~= 0);
   tau1e = tau1e(tau1x)
   tau2e = tau2e(tau2x)
   
   % taus
   initpar(1) = mean(tau1e); % s_lam(3);
   initpar(2) = mean(tau2e); % s_lam(5);
   
   % adjustable parameter initial estimates.
   for i = 1:u(1)
      j = i*3;
      initpar(j) = 5; % offset is min value.
      initpar(j+1) = 100; % 
      initpar(j+2) = 100; % and what remains there
   end;
   
   pmask = [1 1];
   tmax = max(max(x_data));
   % boundaries for values.
   lbound = [initpar(1) initpar(1)];
   ubound = [initpar(2) initpar(2)]; % (not really relevant, since we fix these).
   
   % rest of the parameters.
   for i = 1: u(1)
      ymax = max(y_data(i,:));
      ymin = min(y_data(i,:));
      pmask = [pmask 1 1 1];
      lbound =  [lbound -100 -100 -100];
      ubound =  [ubound  5000 10000 10000];
   end;
   
   order=length(initpar);
   beta = 0;
   alpha = 0;
   model=74;	% see fit_func.m for the implementation
   
   [c,lam]=curve_fitting(x_data(:,10:end), y_data(:,10:end),'simplex','cubic', model, order, initpar,...
      pmask, lbound, ubound, alpha, beta, maxiter );
end;

%now plot it
t = x_data; 
h = findobj('tag', 'multiexp');
if(isempty(h))
   figure('tag', 'multiexp');
else
   figure(h)
end;
clf;

subplot('position', [0.1 0.5 0.8 0.4]);
plot(t', y_data', 'color', 'black');
hold on;
plot(t', fy', 'color', 'green');
if(exist('oldcode'))
   for i = 1:u(1)
      j = i*3;
      z(i,:) = lam(j)+((lam(j+1)*exp(-t(i,:)./lam(1))) + (lam(j+2)*exp(-t(i,:)./lam(2))));   
   end;
end;
plot(t', z', 'color', 'red');


if(exist('oldcode'))
   % return values/arrays.
   tau1 = lam(1);
   tau2 = lam(2);
   offset = lam(3:3:end);
   a1 = lam(4:3:end);
   a2 = lam(5:3:end);
end;

return;



